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STUDEES  IN  STATISTICAL  SIGNAL  PROCESSING 

This  is  a  final  report  on  Contract  AFOSR83-0228,  for  the  period  of  July  1,  1986 
and  June  30,  1988.  Section  I  provides  a  brief  overview  of  our  work,  while  the 
remaining  sections  describe  in  some  detail  our  recent  results  on  efficient  factorization 
of  structured  matrices,  polynomial  zero-location,  handling  of  singularities  in  the 
recursions,  and  recursive  layer  peeling. 

SCIENTinC  PERSONNEL  SUPPORTED  BY  THIS  PROJECT  AND 

DEGREES  AWARDED  DURING  THIS  REPORTING  PERIOD: 

Professor  T.  Kailath 

Graduate  Student:  R.  Ackner 

Senior  Research  Associate:  H.  Lev-Ari 


1.  INTRODUCTION 

The  primary  objective  of  our  research  is  to  develop  efficient  and  numerically 
stable  algorithms  for  nonstationaiy  signal  processing  problems  by  understanding  and 
exploiting  special  structures,  both  deterministic  and  stochastic,  in  the  problems.  We 
also  strive  to  establish  and  broaden  links  with  related  disciplines,  such  as  cascade  filter 
synthesis,  scattering  theory,  numerical  linear  algebra,  and  mathematical  operator  theory 
for  the  purpose  of  cross  fertilization  of  ideas  and  techniques.  These  explorations  have 
led  to  new  results  both  in  estimation  theory  and  in  these  other  fields,  e.g.,  to  new 
algorithms  for  triangular  and  QR  factorization  of  structured  matrices,  new  techniques 
for  root  location  and  stability  testing,  new  recursions  for  orthogonal  polynomials  on 
the  unit  circle  and  the  real  line  as  well  as  on  other  curves,  and  new  approaches  to 
overcome  singularities  and  ill-conditioning  in  the  recursions.  '  ( 

For  several  years,  the  guiding  principle  in  these  studies  has  been  the  concept  of 
generalized  displacement  structure  (Lev-Ari  and  Kailath  (1986)),  which  generalized 
and  subsumed  our  earlier  work  on  Toeplitz-oriented  displacement  structure  (Kailath, 
Kung  and  Morf,  (1979);  see  also  Lev-Ari  and  Kailath  (1984)).  The  same  notion  of 
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displacement  structure  has  also  emerged  from  the  work  of  Heinig  and  Rost 
(1984,1987)  of  East  Germany,  whose  approach  and  methodology  are  significantly 
different  from  ours.  In  particular,  they  focus  only  on  the  problem  of  inversion  of 
structured  matrices  via  algebraic  methods,  while  our  work  has  primarily  addressed 
triangular  factorization  of  such  matrices,  and  our  approach  is  based  on  a  generating 
function  characterization  of  matrices.  In  fact,  in  the  recent  Ph.D.  research  of  J.Chim 
we  were  able  to  reduce  the  inversion  problem  for  structured  matrices  to  the 
factorization  of  certain  block-matrices  with  stmctured  blocks.  This  result  confirms  our 
earlier  findings  (in  Lev- Ari(  1983),  Lev-Ari  and  Kailath  (1984))  on  the  relation  between 
efficient  inversion  and  efficient  factorization  of  structured  matrices:  only  some  of  the 
structured  matrices  that  admit  an  efficient  factorization  procedure  can  also  be 
efficiently  inverted. 

The  generating  function  approach,  which  was  introduced  in  the  Ph.D.  research  of 
Lev-Ari  (1983),  also  suggests  a  natural  geometric  interpretation  of  the  theory,  which 
allows  a  number  of  other  interesting  developments,  e.g.,  studies  of  various  problems  in 
system  theory,  such  as  minimal  realization,  Padd  approximation,  control  design,  and  a 
variety  of  root  distribution  (stability)  problems  for  polynomials.  This  approach  also 
unveils  the  great  flexibility  in  the  computational  details  of  the  factorization  procedure 
for  structured  matrices.  In  contrast,  the  approach  taken  by  Heinig  and  Rost  (and  related 
work  by  Lerer,  Tismenetsky,  Shalom  and  several  other  mathematicians),  leads  to  a 
single  procedure  for  (the  inversion  of)  every  panicular  type  of  structured  matrices.  In 
particular,  our  generating  function  approach  led  us  to  recognize  some  classical  root 
location  procedures,  such  as  the  Schur-Cohn  and  the  Routh-Hurwitz  algorithms,  as 
particular  instances  of  our  factorization  procedure.  Moreover,  by  exploiting  the 
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flexibility  in  our  prototype  procedure  we  obtained  new  alternatives  to  these  classical 
algorithms,  with  reduced  computational  requirements. 

2.  FACTORIZATION  OF  STRUCTURED  MATRICES 

Our  early  work  on  factorization  and  inversion  of  Toeplitz  and  close-to-Toeplitz 
matrices  led  us  to  the  observation  that  for  certain  matrices  the  displacement  matrix 

VjR  :=  R  -  ZRZ'  .  Z  =  [  5,^  ],% 

has  low  rank.  In  particular,  the  displacement  rank  (i.e.,  the  rank  of  V^R)  is  2  for 
both  Toeplitz  matrices  and  for  their  inverses.  We  have  shown  in  previous  work 
(largely  supported  by  AFOSR)  that  the  displacement  concept  is  a  key  tool  for,  . 
developing  fast  algorithms  of  many  kinds,  including  factorization  and  inversion  of 
Toeplitz  and  near-Toeplitz  matrices,  as  well  as  fast  (generalized  Levinson  and  Schur) 
algorithms  for  solving  linear  systems  with  such  coefficient  matrices.  Not  surprisingly, 
these  results  led  naturally  to  cascade  orthogonal  structures  (generalized  lattice  filters) 
for  the  prediction  of  nonstationaiy  processes  (Lev-Ari  and  Kailath  (1984)).  We  have 
also  found  that  the  same  concept  is  tightly  connected  to  the  more  general  problem  of 
cascade  filter  synthesis  in  network  theory  and  digital  filtering  as  well  as  to  a  variety  of 
inverse  scattering  problems  (some  references  are  Rao  and  Kailath,  (1984,  1985), 
Bruckstein  and  Kailath  (1986)). 

Later  we  extended  the  concept  of  displacement  structure  to  a  very  broad  family  of 
structured  matrices,  including  Hankel  matrices  and  their  inverses,  sums  of  Toeplitz  and 
Hankel  matrices  and  several  others  (Lev-Ari  and  Kailath,  (1986)).  The  generalized 
displacement  of  a  matrix  R,  is  defined  as  d(Z,Z)R  where 
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d(A,B)R  :=  Z  ^/t./Z*R(Z*y  ,  (la) 

kj=0 

and  the  asterisk  (*)  denotes  Hermitian  transpose  (complex  conjugate  for  scalars). 
The  generalized  displacement  operator  d(Z,Z)  is  characterized  by  a  (Hermitian) 
displacement  matrix  , 


:={dkj  :  Q<kJ  ^N}  . 


(lb) 


The  previous  notion  of  displacement  corresponds  to  the  particular  displacement  matrix 


Ja  = 


1  0 
.0  -1. 


Jj  , 


while  the  displacement  notion  used  by  Heinig  and  Rost  (1984)  for  Hankel  matrices,'  ' 
corresponds  to  the  displacement  matrix 


We  have  shown  that  efficient  factorization  procedures  for  a  Hermitian  matrix  R  can 
be  formulated  if,  and  only  if,  d  (Z,Z)R  has  low  rank,  and  if  the  displacement  matrix 
has  inertia  (1,1).  This  means  that  must  have  the  form 


h  =  W  -  W*  (2) 

where  <|>,y  are  column  vectors  of  arbitrary  length. 

The  concept  of  displacement  structure  and  its  properties  are  more  conveniently 
described  in  terms  of  generating  functions.  The  generating  function  of  a  matrix  R  is 
a  power  series  in  two  complex  variables,  viz., 

/?(z,w):=[l  z  z^  ...]R[1  w  ...]*  (3) 
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The  displacement  d(Z,Z)R  of  a  matrix  has  the  generating  function  d{z,w)R(z,w), 
where  d(z,w)  is  the  generating  function  of  the  Hermitian  matrix  ,  viz., 

d(z,w)=  X  4,/ (4) 

kj^ 

Thus  the  generating  function  of  a  Hermitian  matrix  with  a  displacement  structure  has 
the  form 


R(z,w) 


G(z)J  GV) 
diz,w) 


(5) 


where  J  is  any  constant  nonsingular  Hermitian  matrix.  The  triple  {d{z ,w),G{z)J  ] 
is  called  a  generator  of  R  (z  ,w ). 

r 

We  have  extended  our  previous  work  (Lev-Ari  and  Kailath  (1984))  on  efficient 
factorization  of  matrices  with  displacement  structure  to  accommodate  the  generalized 
displacement  d(Z,Z)R,  and  we  have  shown  (Lev-Ari  and  Kailath  (1986),  Lev-Ari, 
Bistritz  and  Kailath  (1987))  that  efficient  factorization  of  R  is  possible  if  there  exist 
matrix  functions  0,  (z)  and  constant  matrices  /,•  (Jq  :=J)  that  satisfy  the  matrix 
equation 


0i(z)  A+i  0*(vv) 


Ji- 


d(z,w) 

d(z,C.)d(C/.H') 


JiMiJi 


(6a) 


for  arbitrary  ,  where 


M,  :=  C,‘(Ci)Rr'(Ci.C,)Oi(;i)  =  Mi’  ,  (6b) 

We  have  also  shown  that  (6)  has  a  solution  if,  and  only  if,  /,  are  all  congruent  to  J 
and 


d(z,w)  =  (|)(z)<I>*(w)  -  v(z)y*(w) 


(7) 
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which  is  the  generating  function  version  of  (2).  The  standard  choice  Ci  =  0  produces 
triangular  factorizations;  other  choices  are  required  in  root-location  and  filter  synthesis 
procedures  (see,  e.g.,  Deprettere  and  Dewilde  (1980),  Vaidyanathan  and  Mitra  (1984)). 
The  (nonunique)  solution  of  (6)  is 


@iiz)=il  - 


d(2,Xi) 


where  C/,  is  any  constant  matrix  such  that 

,  (8b) 

and  T,-  is  any  complex  constant  such  that 

d(t,-,T,)  =  0  .  (8c) 

The  factorization  of  R  is  obtained  via  the  recursion 

(z-Ci)G.„i(z)  =  G,(z)0i(z)  i  =  0,1,2,  ...  (9) 

where  0,(z)  have  the  form  (8).  This  algorithm  requires  0(,n^)  computations  to 
factor  a  stractured  n  x  n  matrix,  in  contrast  to  the  conventional  LDL*  algorithm 
which  requires  G(n^)  operations  to  factor  an  arbitrary  n  x  n  matrix. 

It  should  be  observed  that  the  formulation  (8)  allows  a  great  flexibility  in  the 
selection  of  the  scalars  C;  .  'C/  the  matrices  G,  .  The  “universe  of  choice”  has 
four  distinct  dimensions  that  determine  the  specific  form  of  an  efficient  factorization 
procedure  for  a  given  displacement  function  d{z,w)  ; 

(i)  Generator  Type:  If  J  =  Jj  we  say  that  the  generator  {G(z),y}  and  the 
corresponding  recursion  are  of  the  scattering-type.  This  terminology  is  motivated 
by  the  observation  that  with  7  =  we  can  rewrite  (5)  as 


(10a) 


^(z,w)  = 
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u(z)u*(w)  -  v(z)v*(w) 
diz  ,w) 


and  that 


R  >  0 


sup 


V(2) 


z  €  Q»  w  (z  ) 


<  1 


(10b) 


This  means  that  the  ratio  v(z)lu(z)  can  be  interpreted  as  a  (generalized) 
scattering  function  of  a  passive  system.  In  particular,  this  ratio  is  a  continuous¬ 
time  scattering-function  when  d(z,w)  =  z  +  w*  and  a  discrete-time  scattering 
function  when  d(z,w)=  1  -zw*.  Similarly,  if  J  =Jj  :=  gj  we  say  that 

the  generator  {G(z)  ,7}  and  the  corresponding  recursion  are  of  the  immittance- 
type.  In  this  case  (5)  can  be  rewritten  as 


R(z,w) 


which  implies  that 


f(z)g*(w)  +  giz)f  *(z) 
d(z.w) 


(11a) 


R  >  0  w  inf  {Re  4^)  ^  0  . 
zeh.  ^  /(z)^ 


(11b) 


Thus  the  ratio  giz)/f(z)  is  positive  real  m  Q+ ,  which  for  d{z,w)  =  \  -  zw* 
(respectively  z  +  w*)  corresponds  to  a  discrete-time  (respectively  continuous¬ 
time)  impedance/admittance  (=  immittance)  function.  Other  choices  of  J 
determine  other  types. 

(ii)  Extraction  Point:  The  extraction  points  (C/)  can  be  either  on  the  curve  Qq  > 
defined  by  d(z,z)  =  0,  or  in  one  of  the  (disjoint)  domains  Q+ ,  Q_.  In  most 
applications  there  is  an  obvious  preference  for  fixed,  real-valued  extraction  points, 
i.e.,  ^  for  all  / ,  where  Is  real.  However,  in  other  cases  it  may  be 

necessary  to  use  a  different  (complex-valued)  at  each  extraction  step  (see. 
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e.g.,  Vaidyanathan  and  Mitra  (1984)).  Finally,  the  choice  of  affects  the 
computational  complexity  of  evaluating  GjC^,);  the  best  choices  from  this  point 
of  view  are  =  0  and  (probably  next  best)  =  ±1. 

(iii)  Recursion  Depth;  The  fundamental  recursion  (9)  is  a  two-term  vector  recursion. 
This  means  that  each  step  of  the  recursion  relates  variables  indexed  by  i  and  by 
i  +  1  (hence  "two-term")  and  that  there  are  several  functions,  combined  into  a 
vector  G,  (z)  for  each  index  point  i.  In  other  words,  we  can  interpret  (9)  as  a 
linear  vector  difference  equation  of  the  first  order.  Such  equations  can  always  be 
transformed  into  a  single  scalar  linear  difference  equation  of  a  higher  order.  In 
particular,  if  G(z)  has  two  elements  (i.e.,  if  J  is  a  2x2  matrix),  then  (9),  . 
can  be  transformed  into  a  scalar  linear  difference  equation  of  the  second  order, 
also  known  as  a  three-term  recursion.  For  instance,  the  Schur-Cohn  is  a  two- 
term  recursion  for  a  vector  function  of  length  two,  while  the  Bistritz  test  (see 
Bistritz  (1984))  is  a  three-term  recursion  for  a  scalar  function. 

(iv)  Recursion  Parameters:  The  scalars  {x, )  and  the  matrices  {(/,  }  can  be 
chosen  arbitrarily,  subject  to  the  constraints  (8  b),  (8c).  It  seems,  however,  that 
minimization  of  computational  complexity  dictates  specific  (and,  essentially, 
unique)  choices  for  these  parameters. 

Our  recent  research  has  identified  several  choices  that  result  in  recursions  whose 
computational  requirements  are  less  than  those  of  known  procedures. 

In  particular,  we  have  constructed  computationally-improved  alternatives  to  the 
classical  Schur  (and  Levinson)  algorithm,  and  to  the  Schur-Cohn  test,  as  follows: 
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(i)  Immif*  nce-domain  three-term  equivalents  of  the  Levinson  and  Schur  recursions 
for  quasi-Toeplitz  matrices  with  complex-valued  elements  (Bistritz,  Lev-Ari  and 
Kailath  (1987)).  This  extends  our  previous  results  for  real- valued  quasi-Toeplitz 
matrices  (Bistritz,  Lev-Ari  and  Kailath  (1986)). 

(ii)  An  alternative  to  the  Schur-Cohn  algorithm  for  root-location  with  respect  to  the 
unit  circle;  the  form  of  this  new  recursion  is  an  exact  equivalent  of  the  Routh- 
Hurwitz  algorithm  (Lev-Ari  (1987)). 


On  the  theoretical  level  we  have  recently  established  the  following  fundamental 
result:  for  every  d(z,w)  of  the  form  (7)  and  for  every  (finite)  matrix  R, 

In  {d(Z,Z)R)  =  In  {d(Z,Z)[rcv  R-^]}  (12) 


where  the  reversed  matrix  rev  R“^  is  obtained  by  transposing  R~^  with  respect  to 
the  secondary  diagonal.  Thus  the  reversal  of  any  matrix  Q  is  formally  defined  as 


rev  Q  ;=  I  I  (13a) 

where  the  superscript  T  denotes  the  conventional  (non-Hcrmitian)  transpose,  and 


I 


0 


0 


(13b) 


The  fundamental  result  (12)  implies  that  R  and  rev  R~*  have  the  same 
displacement  structure.  It  does  not  tell  us,  however,  how  to  relate  the  parametrizations 
of  generators  of  R  and  of  rev  R“^  Such  relations  are  known  for  Toeplitz  matrices 
and  give  rise  to  the  Levinson  algorithm  and  the  Gohberg-Semencul  formula.  They 
have  been  extended  by  Lev-Ari  and  Kailath  (1984)  to  elose-to-Toeplitz  matrices,  i.e.. 
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to  matrices  with  displacement  structure  involving  d(z,w)  =  1  -  zw*. 

The  proof  of  the  result  (12)  is  based  on  several  observations,  including  the 
following: 

(i)  Under  the  constraint  (7),  the  generalized  displacement  d  (Z,Z)R  is  given  by 

d(Z,Z)R  =  (t)(Z)  R  <1>*(Z)  -  \j;(Z)  R  \i/*(Z) 

(ii)  We  can  assume  without  loss  of  generality  that  ({)(:)  and  \|/(z)  do  not  have 
common  zeros;  if  they  had,  we  could  simplify  d(z,w)  by  extracting  those 
common  zeros.  In  particular,  we  may  assume  that  ())(0)  *  0. 

(iii)  As  a  consequence  of  (i)  and  (ii)  we  conclude  that 

In  {d(Z,Z)R}  =  In  {R-FRF*} 

where  F  :=  (t)"nZ)V(Z),  is  a  lower-triangular  Toeplitz  matrix  whose  first  column 
is  determined  by  the  coefficients  of  the  power-series  expansion  of  \j/(z  )/(})(z ). 

We  have  recently  shown  that  our  factorization  procedure  can  be  extended  to  matrices 
for  which  the  generalized  displacement  R  -  FRF*  has  low  rank,  where  F  can  be 
any  lower  triangular  matrix;  when  F  is  a  lower-triangular  Toeplitz  matrix  the  new 
procedure  reduces  to  the  one  reported  in  [Lev-Ari  and  Kailaih  (1986)].  The  efficiency 
of  this  new  factorization  procedure  depends  upon  the  structure  of  the  matrix  F, 
because  the  procedure  involves  repeated  calculation  of  products  between  certain 
vectors  and  matrices  of  the  form  (F -/,  <!)(! -/,*  F)~‘,  where  /j,,  denotes  the 
(/,/)  element  of  the  matrix  F.  The  computational  complexity  of  such  products 
depends  upon  the  structure  of  F. 


-  12- 


Very  recently,  Mr.  J.  Chun  has  shown  how  F  matrices  composed  Z  of 
different  sizes  can  be  nicely  employed  to  obtain  efficient  algorithms  for  the  least- 
squares  solution  of  various  special  systems  including,  for  example,  triangular  and 
orthogonal  factorization  of  composite  matrices  such  as  (T (7 (T^T)~^T^ . 


3.  BEZOUTIANS  AND  EFFICIENT  ZERO-LOCATION 

One  fascinating  application  of  the  notion  of  generalized  displacement  structure  is 
the  construction  of  Bezoutian  matrices.  Originally,  such  matrices  were  associated  with 
stability  (and  zero-location)  tests  with  respect  to  the  unit  circle  and  the  imaginary  axis. 
We  have  extended  this  notion  to  structured  matrices  whose  generating  function  has  the. 
form 

b(2.w)  =  (14,) 

fl  Zl/4(1  H.]' 

where 


Jd 


(14b) 


and  the  sharp  (#)  denotes  a  suitably  defined  polynomial  transformation.  The 
remarkable  property  of  such  matrices  is  that  most  of  their  elements  vanish  except  the 
elements  in  the  leading  n  xn  principal  submatrix,  where  n  :=  deg  p  (z ).  This 
submatrix,  which  we  denote  by  B,  has  full  rank  (i.e.,  rank  B  =  n)  when  p(z)  and 
p*(z)  are  coprime.  Furthermore,  the  inertia  of  the  Bezoutian  matrix  B  (i.e.,  the 
number  of  its  positive,  zero  and  negative  eigenvalues)  serves  to  locate  the  zeros  of  the 
polynomial  p{z)  with  respect  to  the  following  partition  of  the  complex  plane. 


Q4.  =  {z  ;  d{z,z)  >  0) 
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Cio  ;=  {z  ;  d(z,z)  =  0}  (15) 

:=  {z  ;  d(z,z)  <  0} 

More  specifically  the  inertia  of  B  indicates  how  many  zeros  are  shared  by  p(z)  and 
p*iz)  and  how  many  of  the  remaining  zeros  are  in  £2+  and  in  .  Our  fast 
factorization  procedure  makes  it  possible  to  determine  the  inertia  of  a  Bezoutian  matrix 
in  0(n^)  operations.  For  Bezoutians  on  the  imaginary  axis  and  the  unit  circle  our 
formulation  leads  (among  other  possibilities)  to  the  Routh-Hurwitz  and  Schur-Cohn 
tests,  and  serves  to  delimit  the  family  of  0(/i^)  polynomial  zero-location  procedures. 
Indeed,  we  conjecture  that  every  0(n^)  procedure  corresponds  to  some  specific 
choice  in  the  four-dimensional  “universe  of  choice”  that  we  described  in  Section  2.  . 

4.  SINGULAR  RECURSIONS 

The  fundamental  recursion  step  (9)  can  be  carried  out  if,  and  only  if,  the  scalar 
<10®*  not  vanish.  Thus,  the  recursion  can  be  completed  (for  /  =0,l,...,n) 
if,  and  only  if,  the  matrix  R  is  strongly  regular,  i.e.,  all  its  leading  principal  minors 
do  not  vanish.  Various  solutions  have  been  proposed,  to  overcome  recursion 
singularities  (i.e.,  the  case  when  (^^  )  =  0),  mostly  involving  some  perturbation  of 

the  data  or  the  recursion  parameters.  For  instance,  Vaidyanathan  and  Mitra  (1987) 
suggest  moving  the  point  Ci  >  so  that  i?,  (^,  ,C,  )  ^  0  and  the  recursion  can  be  carried 
on.  This  approach,  in  addition  to  being  rather  ad  hoc,  cannot  be  used  in  applications 
that  involve  predetermined  values  for  {^,},  such  as  the  commonly  used  choice 
Ci  =  0.  There  is  a  variety  of  similar  approaches  in  the  literature. 

Our  approach  to  overcome  recursion  singularities  is  different  and  is  based  on  the 
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observation  that  a  Hermitian  matrix  with  singular  leading  principal  minors  has  a 
block-triangular  factorization  of  the  form  R  =  L  D  L*  where  L  is  lower  triangular 
with  unity  diagonal  elements  and  D  is  a  block-diagonal  matrix.  The  sizes  of  the 
blocks  in  D  are  determined  by  the  rank-profile  of  R,  i.e.,  by  the  variation  of  the 
nullity  (=dimension  of  null  space)  of  the  leading  principal  submatrix  Rq.^  as  a 
function  of  n.  This  implies  that  whenever  /?,  (Ci.Ci)  =  0  it  may  still  be  possible  to 
"jump  over"  the  singularity  and  compute  G,+v(z)  directly  from  G,(z)  via  a 
modified  version  of  (9),  where  v  is  the  size  of  the  block  in  D  whose  top  left 
element  occupies  the  (/  ,i )  position  in  D. 

We  have  recently  derived  this  modified  version  of  the  fundamental  reciu^ion  (9).  . 
for  two  particular  examples:  the  Routh-Hurwitz  algorithm  and  the  Schur-Cohn 
algorithm.  The  former  involves  a  generator  of  the  form 


{d(z,w).  G(z),  J]  =< 


♦ 

z  +  w 


{p(z)+p*i-z)  p(z)-p*(,-z)] 


fo  l] 

U  oj 

J 


where  p{z)  is  a  given  polynomial  (the  recursion  parameters  are 

Ci  =  0,  Xj  =  yoo,  Ui  =  /).  The  latter  involves  a  generator  of  the  form 


{d(z,»v),  G(z),  j) 


1 -zw*  ,  [p(z)+p*(z)  p(z) -p*(z)]  ,  [j  q] 


where  /7*(z)  :=  z***^  [p(l/z*)]  (and  the  recursion  parameters  are 

Ci  =0,  X,-  =  1,  Ui  =J).  Our  modified  recursion  carries  out  the  transformation 
G,  (z)  ->  Gi+^(z)  as  a  sequence  of  single-step  tran^ormations,  viz., 


Gi(z)  — >  G,+i(z)  — >  G,+2(^)  — >  •  •  •  ^  G,+v(z) 


where  the  intermediate  generators  Gi^^iz)  (,i  ^k  -  1),  serve  to  determine  the 
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triangular  factor  L  of  R  column  by  column,  rather  than  by  block-columns  (Pal  and 
Kailath,  (1988)). 

5.  RECURSIVE  LAYER  PEELING 

The  fundamental  factorization  procedure  for  structured  matrices,  viz., 

-  ^i)GM(z)  =  Gi(z)&i(z) 

is,  at  the  same  time,  also  a  layer-peeling,  procedure.  Starting  with  Gq{z),  which  we 
can  interpret  as  boundary  data  for  a  layered  medium,  we  identify  an  elementary  layer 
with  chain-scattering  matrix  0o(z),  then  “peel”  it  off  to  obtain  Gi(z),  the 
boundary  data  for  the  rest  of  the  medium  (with  the  first  layer  removed),  and  repeat  the 
same  procedure  again  and  again.  Such  layer-peeling  recursions  have  been  used  in 
cascade  filter  synthesis  (see,  e.g.,  Dewilde,  Vieira  and  Kailath,  (1978);  Vaidyanathan 
and  Mitra  (1984)),  in  inverse  scattering  (Bruckstein  and  Kailath  (1987)),  zero-location 
(Lev-Ari,  Bistritz  and  Kailath  (1987)),  and  nxxlel-order  reduction  (see,  e.g.,  Genin  and 
Kung  (1981)). 

The  classical  work  of  Schur  (1917)  forms  the  basis  for  much  of  the  subsequent 
work  on  layer  peeling  procedures.  Schur’s  algorithm  associates  a  sequence  of  so- 
caUed  reflection  coeflicients,  all  with  magnitude  bounded  by  unity,  with  every  passive 
scattering  function,  i.e.,  a  function  S(z)  that  is  analytic  and  bounded  by  1  in  the 
unit  disc.  In  particular,  if  S(z)  is  an  all-pass  function,  which  means  that 
1 5  (e-'®)  1=1  for  all  0,  then  Schur’s  algorithm  produces  a  finite  sequence  of 
reflection  coefficients  {/:,■;  O^i^n],  where  |  =1  and  |1:,  |  <1  for 
0^1  ^  n  -  1.  Another  property  of  the  algorithm  is  that  starting  with  a  passive 
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scattering  function  it  generates  a  sequence  of  such  functions.  This  is  the  essence  of 
layer  peeling:  a  single  step  of  the  Schur  algorithm  applied  to  a  passive  medium  leaves 
a  medium  with  the  same  property,  which  makes  it  possible  to  apply  the  same  step 
again  and  again. 

In  addition  to  the  recursive  characterization  of  passivity  via  the  constraint  on  the 
magnitude  of  the  reflection  coefficients,  Schur  also  introduced  an  operator-norm 
characterization  of  passivity:  he  proved  that  for  any  function  /(z)  that  is  analytic  in 
the  unit  disc  we  can  construct  an  infinite  lower-triangular  Toeplitz  matrix  whose  first 

column  consists  of  the  coefficients  of  the  power  series  expansion  of  /  (z ),  viz., 

•  ^ 

/o 

fx  fo  0 

W)=f2fi'  (16a) 

such  that 

sup  |/(z)|^l  <=>  ||L(/-)||2^1  (16b) 

Izl  <  1 

where  ||i4  II2  denotes  the  conventional  (spectral)  norm  of  a  matrix  A ,  i.e., 

IUIl2-=-y“p  |.  II  (16c) 

*  II-*  II2 

and  II X  II 2  denotes  the  Euclidean  (1 2)  norm  of  a  vector  x. 

Schur’ s  algorithm  has  been  later  applied  also  to  functions  with  poles  in  the  unit 
disc,  but  only  to  rational  “all-pass”  functions,  i.e.,  to  functions  /(z)  of  the  form 

/<.)=X^  (17a) 

where  p*(z)  denotes  the  conjugate  reversal  operation,  viz.. 
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p*(z)  =  z‘^«^<*)|>(l/z*)]*  (17b) 

The  well  known  Schur-Cohn  test  associates  with  each  such  functicn^  a  finite  sequence 
of  reflection  coefficients,  some  of  which  have  magnitudes  larger  than  1.  Moreover,  it 
follows  from  the  properties  of  Bezoutians  on  the  unit  disc  that  the  number  of  poles  of 
f  (z)  =  p* iz)tp(z)  inside  the  unit-disc  equals  the  number  of  singular  values  of  the 
matrix  h(f)  that  are  larger  than  |X|  or,  equivalently,  the  number  of  negative 
eigenvalues  of  the  following  finite  rank  matrix, 

R:=  |XpI-L(/ )L‘(/ )  .  (18) 

It  also  follows  that  the  so-called  Pick  matrix 


-  ,  -/(?,)/'«/)  _ , 

R  1 - : — ^  :  0  £  I  j  S  n  1 

1-CiC; 

is  congruent  to  R  of  (18),  viz.. 


ICoa--- 

1  Co  Co'  •  ■  1” 

1  Cl  cf  •  •  ■ 

1  Cl  Cf  •  ■  ■ 

• 

H- 

• 

1  c,  ci  ■  ■  ■ 

^  J 

1  c,  c,'  ■  ■  • 

>  « 

(19) 


and  therefore  has  the  same  number  of  negative  eigenvalues.  The  significance  of  these 
results  has  recently  been  recognized  in  the  context  of  approximation  and  model-order 
reduction  (see,  e.g.,  Genin  and  Kung  (1981)). 

Rational  allpass  functions  of  a  given  degree  k  are  members  in  the  family 
which  consists  of  all  functions  with  k  poles  (or  less)  inside  the  unit  circle,  and  whose 
magnitude  is  bounded  on  the  unit  circle,  i.e.,  sup  l/(z) |  <  <».  It  turns  out  that  the 


^  Atiuming  />  (z )  has  no  zeros  at  z  =  0,  and  applying  (he  algorilhm  to  /  (z  )/X. 
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Schur-Cohn  algorithm  does  not  map  the  family  into  itself,  except  when  ;fe=0. 
This  means  that  this  algorithm  does  not  admit  the  same  layer-peeling  interpretation  as 
the  standard  Schur  algorithm.  Nevertheless,  we  have  found  that  it  is  possible  to 
modify  the  Schur  algorithm  in  such  a  way  that  the  resulting  recursion  indeed  maps  the 
family  into  itself  and,  therefore,  admits  the  same  layer-peeling  interpretation  as 
the  classical  Schur  algorithm.  In  addition,  it  produces  a  transmission-line  model  that  is 
similar  to  the  one  associated  with  the  classical  algorithm. 

Furthermore,  our  modified  recursion  applies  to  every  function  and 

not  just  to  allpass  functions.  Each  layer  in  the  resulting  transmission  line  model  has  an 
indicator  or  ‘sign’.  While  a  positive  layer  maps  into  itself,  a  negative  layer  maps,  . 

into  namely,  it  reduces  by  one  the  number  of  poles  within  the  unit  disc. 

Therefore,  the  number  of  negative  layers  in  the  transmission-line  model  that  is 
generated  by  our  modified  algorithm  equals  the  number  of  poles  that  the  function 
/(z)  has  within  the  unit  disc.  This  is  the  key  idea  in  the  construction  of  efficient 
procedures  for  zero-location  and  model-order  reduction. 

This  result  should  have  some  implications  for  control  systems,  where  the 
underlying  systems  (linearized  about  an  operating  point)  are  often  unstable.  We  are 
exploring  these  implications  and  attempting  to  relate  them  to  results  obtained  in  the 
last  few  years  in  the  active  field  called  //^.-control. 
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